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1. Introduction 

Clear ex am ples of the difficulties associated with applying CFD techniques to ap- 
parently simple problems are provided by computational aeroacoustics and compu- 
tational electromagnetics. These applications require accurate wave propagation 
over long distances for a wide range of frequencies, placing a severe demand on 
numerical algori thms , and raising issues related to efficiency, accuracy, compatible 
space and time treatments, high frequency data, propagation along character- 
istic surfaces, isotropy, stable and accurate artificial boundary treatments, and 
nonrestrictive stability bounds. This paper briefly presents two methods for the 
development of finite difference algorithms which are intended to address these 
issues. High order single step explicit algorithms are possible, and examples with 
up to eleventh order accuracy will be shown. High resolution algorithms in the 
sense of [9] are also possible, with amplification factor and relative phase change 
per tim«» step which are virtually 1 for normal mode frequencies in [0, it] and CFL 
numbers in [0, 1]. If our most accurate algorithm is used to propagate an initial 
periodic sine wave, then after five periods with four grid points per wavelength, 
the TrtATfimnm error is O[10 -6 ], and after five hundred thousand periods with eight 
grid points per wavelength, the maximum error is O[10~ 4 J. High order algorithms 
are relatively more efficient [7], and their relative efficiency tends to increase as 
the error bound decreases, as the simulation time increases, and as the spatial 
dimension increases. Our algorithms show several orders of ma g ni tude difference 
in the number of multiplications required to meet an error bound of O[10~ *] at 
five periods of propagation. Applied computations require boundary conditions. 
A new artifical boundary condition that has been developed by Hagstrom [5] will 
be shown for the linearized Euler equations. This boundary condition is local in 
time and does not require information about the distance or direction of either an 
assumed source or an expanding wave front. The boundary and propagation al- 
gorithms have a consistent derivation and similar properties [4]. Additional issues 
are raised by shocks, but are not addressed here. 


2. Algorithms for the Linearized Euler Equations in ID 

As a first example we will consider the linearized Euler equations in one space 
dimension. The equations for the isentropic case can be written in nondimension- 
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alized form as the system 


du ,,du dp n 

di + M di + fc~ 0, 

dp dp du 

i + < + &=°- 


(i) 


where Af is the constant mean convection Mach number, p is the disturbance 
pressure, and u is the disturbance velocity. System (1) can be diagonalized and 
solved by the method of characteristics to produce the general solution 


u(x, <) =^(ut(x - (Af + 1)<) + p*(x — (Af + l)t)) 
+|(«*(x ~{M - l)t) - pi(x - (Af - l)i)), 
p(x,t) =^(| n'(x - (Af + l)f) + u*(x - (Af + 1 )i)) 
+^(p*(* - (M - 1)<) - ti*(x - (Af - 1)*)), 


( 2 ) 


with intial data u*(x) = u(x,0), and pi(z) = p(x, 0). 

The first step in producing algorithms for solutions of (1) is to locally inter- 
polate u and p at t n with order D polynomials in x: 


D D 

u(Xi + X, t n ) « uo(x) = U «* xa ’ J>( x » + X > *») * P°(*) = P “ xa - 

ar=0 a=0 


The expansion coefficients are obtained by the Method of Undetermined Coeffi- 
cients wring the known data on a given stencil. Note that there is no specification 
of a particular mesh or data type, and that separate derivative terms are not even 
considered. The use of single interpolation polynomials of order D is equivalent 
to simultaneously approximating derivatives up to the D tk order, with 




1 d“ua ^ 1 d a u 
a\~dx S ~ ** a! dx“’ 


and p a = 


1 d a pa _ 1 d a p 
a! dx a ~ a! dx a ' 


( 4 ) 


The local interpolations (3) to u and p at time t n are used as initial data for the 
exact solution (2), producing an approximate local solution 


u(x i + x, t n + 1 ) «i(ua(x - (Af + 1 )t) + pa{x — (Af + 1)<)) 

+i(ua(x - (Af - 1)<) - pa(x - (Af - 1)<)), 

p(xi + x, t n + *) w^(pa(x - (Af + 1)<) + ua(x — (Af + 1)<)) 

+^(pa(x — (Af — l)t) — ua(x — (Af — 1)<)). 
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The approximate local solution (5) is a function of x and t, and it is an exact 
solution to (1) with the local interpolants (3) as initial data, so that (5) is an 
exact local propagator for (1) and correctly incorporates the dynamics of (1). The 
fundamental viewpoint of this method is to approximate the solution of the system 
instead of particular terms in the governing equations. This fundamental idea is 
seen in the use of Riemann solvers for problems with shocks [2], and in the use 
of local Fourier decompositions and separation of variables to develop algorithms 
for a wide variety of problems [1], The approximate local solution (5) is used to 
obtain a computational algorithm at the stencil center, with 

«„ + *)« < +1 = + 1)*) + M-(M + 1)*)) 

+ 1(« a(-(M - l)i) -po(-(M - 1)*)), 

\ («) 

p(x„(. +t)«p" +1 = + 1)*) + ua(-(M + l)t)) 

+ ~ 1 )*) - «<.(-(** - 1 )*)). 

Algorithm (6) uses the exact local propagator (5) with approximate local data (3), 
so that the time evolution introduces no new error, but merely propagates what 
has been introduced by the interpolation. Note that finite difference forms are 
not specified, and that (6) represents a family of algorithms dependant upon the 
interpolant (3). If order D interpolants are used for a particular realization of (6), 
then the algorithm will have accuracy of order D in both space and time. There 
are several interpretations of algorithm (6). It may be viewed as an application of 
the Method of Characteristics, since in this case, the general solution form for the 
exact local propagator is obtaind by this method. If the interpolation expansion 
coefficients are viewed as spatial derivatives (4), then the local solution in space 
and time (5) can be viewed as a Cauchy-Kowaleskaya expansion. The algorithm 
can also be reformulated as a truncated Taylor series expansion in the time step 
size k , with 

u" +1 =u 0 - k(pi + Mu x ) + k 2 (2Mpi + (M 2 + l)u 2 ) 

- Jb 3 ((l + ZM 2 )pz + M( 3 + M 2 )u 3 ) + • • • , 

p? +1 =Po - k( Ul + Mpx ) + k 2 (2Mu 2 + (M 2 + 1 )P2) C 

- Jb 3 ((l + 3 M 2 )u 3 + M( 3 + M 2 )p3) + . . . , 

where the grid ratio A = £ is implicit in (7), since the coefficients u a and p a 
include the factor h~ a for space step size h. A truncated Taylor expansion in time 
can be viewed as a generalized Lax-Wendroff method [8]. Algorithm (6) can also 
be reformulated as a conventional finite difference method, since the underlying 
spatial interpolations use local polynomials. Further details are in [3]. 
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We will introduce four relatively conventional realizations of algorithm (6), 
with central stencils for the interpolants (3), and with values for u and p at each 
grid point. These algorithms have three, five, seven, and nine point central sten- 
cils, they are second, fourth, sixth, and eighth order accurate in both space and 
time, and they are refered to as the ”c3o0ex,” ”c5o0ex,” ”c7o0ex,” and n c9o0ex” 
methods, respectively. These four algorithms are all single step explicit methods 
with dispersive truncation errors, and each is stable for j < y^ jf- Grid refinement 
data for these algorithms is presented in Figure 1, confirming the order of accuracy 
of each of the four methods. This grid refinement data is discussed in more detail 
with data from an additional class of high resolution Hennitian methods that are 
introduced below. 

Stable high order boundary treatments can be derived if these algorithms are 
viewed as C auchy- Kowaleskaya expansions. The interpolant from the data on a 
single stencil next to a boundary is used as initial data for the approximation 
of the evolution of the solution over the interval from the stencil center to the 
boundary point. Values for u and p are computed at each grid point that is not 
on the boundary in this interval, and outgoing Riemann variables are computed 
at the boundary point. The boundary treatment for the c9o0ex method must be 
modified for stability, and it uses a truncated eight point boundary stencil with 
additional data at the boundary point. Table 1 presents results for a simulation 
with inital data t*(x, 0) = 0 and p(x, 0) = Sin(xx), for — 1 < x < 1, with A = 0.8 
and M = 0. For the boundary treatment, ti+p is computed at the x = 1 boundary, 
u — p is computed at the x = — 1 boundary, and typical CFD boundary conditions 
are imposed, with u(— l,t) and p(l,t) given. The data in Thble 1 shows that the 
propagation algorithms and boundary treatments are stable with from second to 
eighth order accuracy in both time and space. Further details are in [3]. 


TABLE 1: Data From Explicit Algorithms with Boundary Treatments 
Maximum Error in u or p at t = 10 


2 

h 

nio 

c3d0ex 

c5d0ex 

c7d0ex 

c9d0ex 

8 

50 

1.87D-01 

1.17D-02 

1.01D-02 

8.65D-05 

16 

100 

4.57D-02 

9.98D-04 

5.26D-05 

8.50D-07 

32 

200 

1.32D-02 

8.11D-05 

8.79D-07 

4.69D-09 

64 

400 

3.50D-03 

5.54D-06 

1.29D-08 

2.09D-11 

128 

800 

8.94D-04 

3.58D-07 

1.89D-10 

3.75D-13 

256 

1600 

2.25D-04 

2.27D-08 

3.35D-12 

4.91D-13 

512 

3200 

5.66D-05 

1.43D-09 

1.89D-12 


1024 

6400 

1.42D-05 

9.12D-11 




3. Hermitian Algorithms 

We will introduce a second family of algorithms for (1), which use the exact lo- 
cal propagator (5), but which are distinguished from the relatively conventional 
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algorithms introduced above by the use of Hermiti&n interpolates for (3). This 
particular family of algorithms uses and computes values for u and p, and for 
their spatial derivatives. Various orders of spatial derivative data are required 
at each grid point, depending upon the particular algorithm. If these algorithms 
are viewed as approximate local solutions for u and p, or as Cauchy- Kowaleskaya 
expansions that are locally defined as functions of x and t, then they can be 
differentiated in x to provide local solutions for the spatial derivatives that are 
consistent with the solutions for u and p. The local derivative solutions obtained 
by differentiating (5) in x are used to obtain computational forms for the spatial 
derivatives that are analogous to and consistent with the computational forms (6) 
for u and p. These forms may be expressed as Taylor series time expansions in k, 
with 

ux”"*" 1 =ttj — 2k(p2 + A/t*2) 3fc 2 (2Afp3 + (Af 2 -I- lju^) 

— 4fc 3 ((l + 3Af 2 )p4 + Af(3 + M 2 )u 4) + . . . , 
px" +1 =pi - 2k(u2 + Mpi) + 3fc 2 (2Aftt3 + (Af 2 + l)p3) 

- 4* 3 ((1 + 3 Af 2 )u 4 + Af(3 + Af 2 )p 4 ) + . . . , 

and with similar forms for second and higher derivatives. The interpolation step 
for the Hermitian algorithms only produces interpolants (3) for the functions u 
and p, using local polynomial approximations and the Method of Undetermined 
Coefficients with all of the data at each grid point on a given stencil. Note that 
there is no separate consideration given to interpolation or propagation for the 
derivative data. 

Hermitian algorithms with central stencils have stability problems if the max- 
imum order of accuracy is obtained. Stable algorithms with exact propagators can 
be obtained on alternating grids offset by a half mesh width. We will introduce 
four algorithms which use staggered two point grids, and three algorithms which 
use staggered four point grids. On staggered two point grids, if just u and p are 
used, then a first order exact propagator algorithm is obtained, which we will de- 
note by ”c3o0s2.” If up to first, second, or third derivatives are used in addition to 
u and p, then stable algorithms are obtained which are third, fifth, or seventh order 
accurate in space and time, referred to as the ”c3ols2,” ”030282,” and ”c3o3s2” 
methods, respectively. On staggered four point grids, if just u and p are used, 
then a third order exact propagator algorithm is obtained, which we will denote 
by ”c5o0s2.” If up to first or second derivatives are used in addition to u and p, 
then stable algorithms are obtained that are seventh or eleventh order accurate in 
space and time, referred to as the ”c5ols2,” and ”c5o2s2” methods, respectively. 
These seven algorithms are all single step explicit methods on staggered grids with 
diffusive truncation errors, and each is stable for j < where the staggered 

half step size is j. Note that two half time steps can be composed to produce 
single step algorithms with time step size k. The full time step algorithms result- 
ing from staggered two point grids can be formulated as single step methods on a 
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central three point stencil, and the algorithms from staggered four point grids as 
single step methods on a central seven point stencil. Other combinations are pos- 
sible, such as composing the c3o3s2 and c5ols2 methods, to produce a single step 
explicit seventh order Hermitian method on a single central five point stencil. It 
will be shown below that some of these Hermitian algorithms have extraordinary 
resolution and accuracy. 

4. Numerical Experiments with the ID Algorithms 

Numerical experiments will be conducted with the eleven algorithms that have 
been introduced above, by computing solutions to (1) with inital data u(x, 0) = 0 
and p(z, 0) = Sin(xz) for —1 < x < 1, and with A = 0.8 and M = 0. Periodic 
boundaries will be imposed, withp(— l,t) = p(l,t) and u(— l,t) = u(l,t) for 0 < t. 
Note that the initial data has wavelength and period 2. Figure 1 shows the logio 
of the maximum absolute error at t = 10 in u or p versus the log? of the number 
of grid points per wavelength. This data corroborates the order of accuracy of 
each of the methods. Note that the first order c3o0s2 method is incapable of 
producing accurate results with any reasonable resolution, and that the second 
order c3o0ex or Lax-Wendroff method requires 64 grid points per wavelength to 
produce O[10~ 2 ] accuracy after five periods. Note also that the c7o0ex and cdoOex 
central methods show a greater sensitivity to roundoff errors than the high order 
Hermitian methods. An interesting feature of Figure 1 is the data at the coarsest 
resolution, with 4 grid points per wavelength. At this resolution, the simulations 
with the c3o2s2, c3o3s2, c5ols2 and c5o2s2 algorithms have errors that range from 
0[1O -2 ] to O[10 _# ], in contrast with errors that range from 0(1] to O[10 l ] for 
the other algorithms. An O[10 -6 ] error after five periods of propagation with 
four grid points per wavelength is exceptionally high resolution. Notice also the 
relative errors from methods which are of similar order. The c5o0s2 and c3ols2 
algorithms are both third order, but the c3ols2 algorithm produces lower errors 
by about one order of magnitude at each grid resolution. The c5ols2 and c3o3s2 
algorithms are both seventh order, but the c3o3s2 algorithm produces lower errors 
by about two orders of magnitude at each grid resolution. The conventional sixth 
order c7o0ex and eighth order c9o0ex methods have larger errors at each grid 
level than the fifth order c3o2s2 and seventh order c5ols2 Hermitian methods, 
respectively. In these comparisons of algorithms with similar orders of accuracy, 
the algorithm which produces lower errors has higher resolution from using more 
derivative information at each grid point. The accuracy of a numerical algorithm 
is determined by both its order of accuracy and its resolving power. 

The periodic problem which produces Figure 1 is also used for propagation 
with 8 grid points per wavelength out to t = 10, i = 1,000, and i = 100,000. 
The data is presented in Table 2, where 0[1] errors are marked by asterisks, and 
where the algorithms are ranked by order. Note in Table 2 for t = 10, that the 
error data from the sixth order c7o0ex method is two orders of magnitude higher 
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than from the fifth order c3o2s2 method, and that the error from the eighth 
order c9o0ex method is two orders of magnitude higher than from the seventh 
order c5ols2 method, and four orders of magnitude higher than from the seventh 
order c3o3s2 method. A similar comparison is also seen in Table 2 for the data 
from t = 1,000 and t = 100,000. These comparisons once again show that both 
the order of accuracy and the resolution of a numerical algorithm determine its 
accuracy. Note that the only algorithms which produce errors that are not 0[1] 
at < = 100, 000 are the high order and high resolution c5ols2, c3o3s2 and c5o2s2 
methods. These results also show that far field may be redefined by several orders 
of magni tude, and that efficient propagation to a truly far field requires methods 
which are both high order and high resolution. 


TABLE 2: Long Time Simulations For Each Algorithm 
Maximum Error in u or p at Various Times 


Method 

Order 

II 

I-* 

o 

t = 1,000 

< = 100,000 

c3o0s2 

1 

******** 

******** 

******** 

c3o0ex 

2 

6.76D-01 

******** 

******** 

c5o0s2 

3 

2.42D-01 

******** 

******** 

c3ols2 

3 

3.97D-02 

******** 

******** 

c5o0ex 

4 

9.14D-02 

6.74D-01 

******** 

c3o2s2 

5 

2.27D-04 

2.26D-02 

******** 

c7o0ex 

6 

1.11D-02 

1.04D-01 

******** 

c5ols2 

7 

4.52D-05 

4.55D-03 

3.75D-01 

c3o3s2 

7 

6.74D-07 

6.78D-05 

6.75D-03 

c9o0ex 

8 

1.39D-03 

1.38D-02 

******** 

c5o2s2 

11 

9.33D-10 

9.37D-08 

9.37D-06 


The superior accuracy of the high order and high resolution methods is ob- 
tained by using more complex algorithms, with more variables and equations, and 
with more operations per grid point per time step. It is natural to ask wether 
or not the simpler, more conventional algorithms are more or less efficient than 
the Hermitian algorithms. We will address this issue by replotting the data from 
Figure 1, this time with the log 10 of the maximum error at t = 10 on the vertical 
axis, and the log lQ of the total number of multiplications required for the entire 
simulation out to t = 10 on the horizontal axis. This data is presented in Figure 
2. It has been shown in [7] that relative computational efficiency increases with 
order of accuracy, and that the relative efficiency increase as the error tolerance is 
lowered, and as the simulation time is increased. The data in Figure 2 from our 
methods clearly shows that for every level of error the efficiency in meeting that 
error tolerance increases with the order of accuracy of the algorithm, in spite of 
the fact that algorithms of radically different type are being used. In particular, 
the Hermitian algorithms not only use more degrees of freedom of data in order to 
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inc rease order of accuracy, but they also solve more equations for more variables, 
and yet they are still more efficient than the less accurate algorithms. 


5. Normal Mode Analysis in ID 

The special properties of the full discretization Hermitian algorithms can be shown 
by a normal mode analysis. We will consider Hermitian algorithms applied to the 
linear first order wave equation 


du 

~dt 


+ M 


du 

dx 


= 0 . 


( 9 ) 


The local exact propagator is defined by the Method of Characteristics. A normal 
mode for (9) can be written in local coordinates as 


u(x,t) = aExp[t0- — J- — -], 


( 10 ) 


with amplitude a, frequency 9 € [0, x], and with space mesh size h. The symbol 
p for the algorithm can be written in terms of the normal mode (10), with 


P( A,*) 



( 11 ) 


where u* +1 is obtained from the algorithm with the normal mode u* as the known 
solution at <», and where A = ^ is the CFL number. Spatial derivative data 
is needed by the Hermitian algorithms at each grid point, and since the normal 
mode (10) is perfectly known as a function of x at we will simply take and 
use its derivative values for this data. This heuristic procedure is intended to 
obtain qualitative insights about the algorithms, and is not intended as a rigourous 
stability analysis, which would treat the total system of u and its spatial derivatives 

with independant error modes expected in each. IVam (11) we obtain the norm 

/ 

|| P (A,»)|| = (fi e W 2 + /mW J ) i . (12) 

and the phase change per time step 

pc(A,9) = Cos (13) 


Note that we have varied from [10] by using a definition of phase change in terms 
of Cos -1 rather than Tan -1 . For the phase speed properties of an algorithm, we 
will use the normalized relative phase change per time step 

( 14 ) 
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rather than the phase change per time step (13). 

All of the Hermitian methods with two point staggered grids were applied to 
the wave equation (9), and single step forms of the algorithms were obtained on 
three point central stencils. The norm of the amplification factor (12) and the 
relative phase change per time step (14) of these single step algorithm forms have 
been obtained numerically from an analytical expression. The most accurate of 
these four methods is the seventh order high resolution c3o3s2 algorithm, which 
uses u and its first through third spatial derivative at every grid point. The norm 
of the amplification factor (12) for the c3o3s2 algorithm is plotted in Figure 3(a) 
as a function of the wave number 0 6 [0, r] and the CFL number A € [0, 1]. Note 
t.hnf. Figure 3(a) shows the norm of the amplification factor as less than or equal 
to 1 in the specified parameter range. The most dissipated behaviour is in the 
limi t at 6 = x, and the norm of the amplification factor at 0 = x is plotted as 
a function of A in Figure 3(b). Figure 3(b) shows the norm of the amplification 
factor varying between approximately 0.9988 and 1 at 0 — x. The relative phase 
change per time step (14) for the c3o3s2 algorithm is plotted in Figure 4(a), also 
as a function of the wave number 0 € [0, x] and the CFL number A 6 [0, 1]. The 
most dispersed behaviour is in the limit at 0 — x, and the relative phase change at 
0 = x is plotted as a function of A in Figure 4(b). Figure 4(b) shows the relative 
pknqc chang e per time step varying between approximately 0.9997 and 1.0002 at 
0 = n. The amplification factor and relative phase change per time step plots in 
Figures 3 and 4 are extraordinary, and show truly "spectral like” qualities for the 
c3o3s2 algorithm. Recall from the numerical experiments reported above, that the 
c 5o2s ? algorithm is of even higher resolution than the c3o3s2 algorithm. 


6. Linearized Euler Equations in 2D 

The algorithms that have been developed in one space dimension have all used 
the Method of Characteristics to obtain an exact solution form that propagates a 
local spatial inerpolant. The Method of Characteristics cannot be used in multiple 
space dimensio ns for nondiagonalizable hyperbolic systems. We will consider the 
Euler equations in two space dimensions in order to indicate an approach 
for developing exact local propagators in multiple space d ime nsi o ns. Consider the 
linearized Euler equations for the nondimensionalized isentropic case, 




(15) 


where (17, V) is the constant mean convection velocity in Mach number, p is the 
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pressure disturbance, and (u, t>) is the velocity disturbance. This two dimensional 
system is nondiagonalizable , with wave propagation along characteristic surfaces. 

We will briefly describe the development of a second order explicit algorithm 
for (15) on a symmetric 3x3 stencil. A second order local interpolant to u in i 
and y at i n can be written as 

u(xi + x , yj + y, t„) « ua(x, y) =u 0 ,o + «i,o* + « 2 ,o* 2 

+(u 0 ,i + *1,1* + u 2 ,ix 2 )y (16) 

+(t* 0 ^ + u M x + u 2 ,2* 2 )y 2 , 

with wwiilw interpolants far v and p. The expansion coefficients are simultaneously 
obtained by the Method of Undetermined Coefficients, and can be interpreted as 
spatial derivatives, with 


1 d° + ^uo 1 d a+ ^u . 

Notice that there are up to fourth order cross derivative terms in the interpolant 
(16). Exact polynomial solution forms for the linearized Euler equations (15) can 
be derived by substituting the expansion forms 



2 4 

u(x< + x, yj + y, t„ + 1) « uo(x, y, t) = ^ 

a, 0=0 T=0 
2 4 

v(x, + X, yj + y, tn + 1) » va(x, y, t) - ^ v a , 0 n x o, y fi t' 1 , (18) 

a, 0=0 7=0 
2 4 

p(x j + x,y i + y,<* + <) «pa(x,y,f)= J ^Paj^ift' 1 , 

a t /J= 0 7=0 

into system (15), and obtaining all the terms with 7 £ 0 by requiring system (15) 
to be satisfied for all x, y and t. Coefficients with 7 ^ 0 are equivalent to time 
derivatives, and the resulting polynomial solutions are expres s ed entirely in terms 
of the spatial expansion coefficients. Note that there will be third and fourth or- 
der time terms present in the exact propagator solution forms. This procedure 
for obtaining an exact solution form is equivalent to the Cauchy-Kowaleskaya pro- 
cedure [6]. The exact polynomial solution forms from (18) give exact propagator 
algorithms with correct multidimensional wave dynamics for the local spatial inter- 
polants (16) as initial data. Since a biquadratic interpolant is used, the resulting 
algorithm will have 0[h 2 + Jk 2 ] truncation error, and it is dispersive. This pro- 
cedure has been used with other symmetric stencils and interpolants to produce 
algorithms with fourth and sixth order accuracy in both space and time (see [3] for 
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details). Because a local exact propagator is used to develop these algorithms, the 
correct propagation of information along characteristic surfaces is automatically 
incorporated in the local approximation of the solution, so that this approach can 
be said to generalize the Method of Characteristics to nondiagonalizable systems 
in multiple space dimensions. 

We have implemented new boundary conditions developed by Hagstrom [5] 
as algorithms that are compatible with our propagation methods for the linearized 
Euler equations in two space dimensions. Hagstrom actually provides a sequence 
of local approximations on the artifical boundary that are defined by a variable 
number of auxiliary functions. Each of the approximate problems are strongly well 
posed, and the approximate solutions converge exponentially to the exact solution 
on the open domain as the number of auxiliary functions increases [5]. As an 
example of these boundary conditions, consider the case of a subsonic mean flow 
in the positive x direction, with normal Riemann variables r = u+p and / = u - p. 
The function r is outgoing at an artifical boundary on the right, and is viewed 
as being defined on the interior and the surface, while the function l is incoming 
and is viewed as being defined only on the surface. The boundary surface values 
of r and t> are obtained essentially in the same way as the outgoing Riemann 
variables in the one dime nsional case considered above, with the algorithm forms 
being interpreted as C auchy- Kowaleskaya expansions in space and time, and all 
relevant solution values obtained over a common boundary stencil. The function 
l is obtained with auxiliary functions fj and gj, from the system 


dt dy 


5^(/i + 9j) = 
>=i 

dfi 

+a M 

Pi ~ 0 

dt 

+ ‘dy 

1 

to 

dgj 

dgj 

ft <?(( - r ) 

dt 

3 dy 

N 

£ 

CS 

II 


(19) 


where ay = y/i — t/ 2 cos( 2t ^* y ) + V , and fy — — sin 2 ( 2 m+i )• Note that sys- 

tem (19) is forced by the interior solution propagating across the artifical boundary 
in the form of the §J term in the equation for /, and the fp terms in the equations 
for the auxiliary functions. Note also that this boundary system does not require 
assumptions about solution form or source location. In practice, we have used 
this condition with disturbance data entirely contained within the computational 
domain, and with boundary data initialized as 0. Similar systems are defined on 
the left hand artifical boundary (see [4], [5] for details). Boundary condition (19) 
and its left hand analog have been implemented in both second and fourth order 
algorithms, and numerical experiments with m = 2 auxiliary functions have shown 
no visible evidence of reflection (see [4] for details). 
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7. Conclusions 

Single step methods with high order accuracy in both space and time are shown 
and compared. A particular class of methods using Hermitian interpolation on 
alternating grids has shown both high order accuracy and spectral like high res- 
olution. The hi gh order and high resolution methods are more efficient than the 
less accurate methods by orders of magnitude, even though the high order and 
high resolution methods are considerably more complex. Calculation to the ‘Tar 
field” ran be redefined as propagation to more than 0(10*] wavelengths or periods. 
Both high order accuracy and high resolution methods are required to compute to 
a true far field, with an example algorithm producing OJlO -4 ] errors after 5 x 10 5 
periods of propagation with eight grid points per wavelength. Algorithms for the 
linearized Euler equations in 2D are discussed, with propagation along characteris- 
tic surfaces, which generalizes the Method of Characteristics to nondiagonalizable 
Hyperbolic systems. Unobtrusive artificial boundary conditions are indicated. 
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Figure 1: Maximum Error by the Number of Grid Points per Wavelength. 



Figure 2: Maximum Error by the Total Number of Multiplications. 
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Figure 3(a): Amplification Factor for the c3o3s2 Algorithm. 


c3o3s2 Algorithm Amplification Factor at theta = Pi 



lambda 


Figure 3(b): Amplification Factor at 6 = r for the c3o3s2 Algorithm. 
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Figure 4(a): Relative Phase Change for the c3o3s2 Algorithm. 



c3o3s2 Algorithm Relative Phase Change at theta = Pi 



Figure 4(b): Relative Phase Change at 6 = ir for the c3o3s2 Algorithm. 
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